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By choosing appropriate deformed Maxwellian ion and electron distribution functions depending 
on the two particle constants of motion, i.e. the energy and toroidal angular momentum, we reduce 
the Vlasov axisymmetric equilibrium problem for quasineutral plasmas to a transcendental Grad- 
Shafranov-like equation. This equation is then solved numerically under the Dirichlet boundary 
condition for an analytically prescribed boundary possessing a lower X-point to construct tokamak 
equilibria with toroidal sheared ion flow and anisotropic pressure. Depending on the deformation 
of the distribution functions these steady states can have toroidal current densities either peaked 
on the magnetic axis or hollow. These two kinds of equilibria may be regarded as a bifurcation in 
connection with symmetry properties of the distribution functions on the magnetic axis. 

PACS numbers: 52.65.Ff, 52.55.-s, 52.55.Fa 


I. INTRODUCTION 

Kinetic equilibria may provide broader and more 
precise information than multifluid or MHD equilibria 
for axisymmetric plasmas governed by the well known 
Grad-Shafranov (GS) equation. Because solving self- 
consistently the kinetic equations is difficult, particularly 
in complicated geometries, the majority of kinetic equi¬ 
librium solutions have been restricted to one-dimensional 
configurations in plane geometry i-i- Of particular in¬ 
terest are equilibria with sheared flows related to electric 
fields which play a role in the transition to improved con¬ 
finement regimes of magnetically confined plasmas as the 
L-H transition and the formation of Internal Transport 
Barriers. Construction of kinetic equilibria is crucially 
related to the particle constants of motion, upon which 
the distribution function depends. In the framework of 
Maxwell-Vlasov theory only two constants of motion are 
known for symmetric two-dimensional equilibria, i.e. the 
energy, E, and the momentum conjugate to the ig- 
norable coordinate X 3 . Therefore for distribution func¬ 
tions of the form f{E,Cx 3 ), only macroscopic flows and 
currents along the direction associated with X 3 can be de¬ 
rived, e.g., toroidal flows for axisymmetric plasmas. The 
creation of poloidal flows requiring additional constant (s) 
of motion remains an open question. 

In the presence of magnetic fields kinetic equilibria con¬ 
structed in the literature usually concern neutral non¬ 
flowing plasmas in connection with the choice of func¬ 
tionally identical electron and ion distribution functions 
in addition to the quasineutrality condition implying van¬ 
ishing electric fields i-B- On physical grounds the as¬ 
sumption is oversimplifying because, in addition to the 
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aforementioned importance of electric fields, it ignores 
the mass difference of ions and electrons. A more realis¬ 
tic treatment permitting finite electric fields consists in 
using only the quasineutrality to express the electrostatic 
potential, $, in terms of the components of the vector po¬ 
tential involved 1-0: e.g. for two-dimensional equilib¬ 
ria in plane geometry, ^{x,y) can be expressed in terms 
of Az{x,y). In connection with the present study we 
also refer to a previous paper [l^ . in which it has been 
proved that the current on the magnetic axis of an ax¬ 
isymmetric Vlasov equilibrium vanishes if the gradient of 
the distribution function and the electric filed are taken 
equal to zero on that axis. However, for the translational 
symmetric two-dimensional configurations, quasineutral 
equilibria with non-vanishin g cu rrent density were explic¬ 
itly found in references |10l|-[lll|. 

The aim of the present contribution is to construct 
Vlasov quasineutral equilibria in toroidal axisymmetric 
geometry. The major new ingredient compared with 
plane equilibria is toroidicity which plays an important 
role in tokamaks. A first step to the construction of such 
equilibria was made in [l3| for deformed Maxwellians 
leading to GS-like equations describing either static equi¬ 
libria or equilibria with rigid toroidal flow. However, 
complete construction of specific equilibrium configura¬ 
tions by solving the GS-like equations was not made in 
[ 1 ^ . Here we select another form of exponentially de¬ 
formed Maxwellians with exponent depending quadrat- 
ically on the toroidal angular momentum. This choice 
amenable to analytic integrations in the velocity space 
results to more pertinent tokamak equilibria with sheared 
toroidal flows. The distribution functions chosen have fi¬ 
nite gradients on the magnetic axis and provision is made 
so that the electric field vanishes thereon. This is impor¬ 
tant because otherwise the resulting E x B drift on axis 
would create unphysical perpendicular flows. Also, de¬ 
pending on the symmetry properties of the distribution 
functions on the magnetic axis determined by appropri¬ 
ate free parameters, derivation of equilibria with toroidal 
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current densities either peaked on axis or hollow is possi¬ 
ble. It is shown that the above described procedure leads 
to a GS-like equation with a transcendental right-hand 
side which is then solved numerically for a fixed diverted 
boundary. Up to the best of the authors knowledge such 
Vlasov tokamak equilibria in toroidal geometry are con¬ 
structed for the first time. 

The organization of the paper is as follows: In Sec. 

II we present the general setting of the equilibrium 
equations in the framework of Vlasov theory. In Sec. 

III making the aforementioned choice of distribution 
functions we derive the GS-like equation. In Sec. IV 
we examine the equilibrium properties by calculating 
macroscopic quantities as the toroidal current density, 
fluid velocities and pressure. Finally the main conclu¬ 
sions are summarized in Sec. V. 


II. GENERAL DEVELOPMENT OF THE 
KINETIC FRAMEWORK 

We consider axisymmetric toroidal plasmas and em¬ 
ploy cylindrical coordinates {z,R,(p), with z correspond¬ 
ing to the axis of symmetry and 4> being the toroidal, 
ignorable angle. The coordinate system is illustrated in 
Figure 6.1 of reference [T3 |. Axisymmetry means that 
any quantity depends solely on z and R. The toroidicity 
relates to the scale factor 1/R appearing in the various 
equations, e.g. equation (IT6)) for the magnetic field; also 
the differential operator in connection with Ampere’s law 
is the elliptic operator on the LHS of Eq. (ITOl) in Sec. Ill 
while it is the Laplace operator in plane geometry. As al¬ 
ready mentioned in Sec. I in axisymmetric plasmas only 
two constants of motion are known, the energy E and the 
angular momentum C. Using convenient units by setting 
the magnetic permeability of vacuum equal to unity, we 
have for the ions 

E, = e<i>{R,z) +YivR + vl + vl) 

Ci = MRv^ + eRA^ ( 1 ) 

while for the electrons we have 

Ee = -e<S>{R,z) + ^{v% + vl+vl) 

Ce = mRv^ — eRA^ ( 2 ) 

where $ is the electrostatic potential, A^ is the toroidal 
component of the vector potential, M and m are the 
masses of ions and electrons respectively, vr, v^, Vz are 
the components of the particle velocities along the basis 
vectors of the cylindrical coordinate system. The charge 
e is taken as the absolute value of the electron charge. 

The solutions of the ion and electron Vlasov equations 
are given as 


( 3 ) 


with a normalization of fi , fe equal to the total number 
of particles N so that the densities are given by 

Ui = J fi{Ei,Ci)(fv 

Ue = J feiEe,Ce)(fv (4) 

The electric current density is given by 

J 4 , = e J v^fi{Ei,Ci)dPv - e J v^fe{Ee,Ce)(fv (5) 

We assume quasineutrality through the whole plasma 
instead of Poisson equation for the electric potential. 
Also we demand that the electric field, E, vanishes on 
axis. Otherwise, as already mentioned in Sec. I, the R 
component of the electric field in combination with the 
purely toroidal magnetic field on axis would lead to an 
unphysical E x B drift parallel to the axis of symmetry. 
A similar condition of vanishing E on axis was adopted in 
[H for a near axis consideration of the Vlasov equation. 
This leads to 


m = He (6) 

everywhere inside the plasma and in particular 

Vrii = Vue (7) 

on the magnetic axis. We introduce now 4/ = RA^ as the 
poloidal flux around the magnetic axis, a function which 
labels the magnetic surfaces and can be taken equal to 
zero on the magnetic axis. Also, we compute explicitly 
both sides of Eqs. and © using (IH) to obtain 


J ME,,a)d^V = J fe{Ee,Ce)d^ 


( 8 ) 


and 


Vm = 


Vrie = 






'1” 


' dfe dfe 


dfe „ p 

V K 


da 


(9) 


We see from Eq. ([5]) that the electrostatic potential will 
be, in general, a function of dt and R in contrast with 
the scale factor free two-dimensional case treated in ref¬ 
erences MM- Similarly Eqs. m show that the gra¬ 
dient of the electrostatic potential does not necessarily 
vanish at the magnetic axis given by 4' = V4) = 0 

because and cannot vanish unless the toroidal 
current vanishes also. This leads us to look for special 
choices for fi, /e which would allow us to have in gen¬ 
eral different from zero, while at the same time satisfy 
dt = Vdt = V$ = 0 on axis. 


f^ = ME,, a) 

fe = fe{Ee.Ce) 
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With this in mind we observe that in order to satisfy 
Eq. ([7]) at the magnetic axis, a sufficient condition would 
be that the equality of the last terms of Eqs. is sat¬ 
isfied at d/ = 0. We thus choose 

fi = no expi-PiEi) x 

1 -f tticxp 


fe = no 


mPe 

2t: 


3/2 


eXp{-PeEe) X 




1 -f UeCXp ( Pl-^{Ce + Bf 


( 10 ) 


with Ps = l/pkTs), (s = i,e) where E^, Ci, Eg, Cg 
are given by Eqs. o, @ and Q < ai, Ug < 1 are 
constants. Also Wc/, Rn are constants to be deter¬ 
mined below. The free parameters A and B are crucial in 
determining the equilibrium characteristics. Specifically, 
for A = B = 0 the distribution functions become sym¬ 
metric with respect to on the magnetic axis (where 
dt = 0) resulting in vanishing toroidal current densities 
thereon and therefore hollow current density profiles. In 
contrast, peaked current density profiles are derived for 
A 0 or/and B 0, a. choice which breaks the v^- sym¬ 
metry of fi and fg on axis. To derive peaked on axis 
and vanishing on the boundary as well as macroscopic 
ion flow we have chosen 

Ap,^ :=A=-1 and Bpg^ := A =-I (11) 

Irrespective of the values of A and B the integrations in 
the velocity space can be performed analytically. The 
characteristics of both kinds of equilibria, peaked and 
hollow, will be presented in Sec. IV. However, since for 
A ^ 0 and B 0 certain expressions become rather com¬ 
plicated, in the reminder of the report complete analytic 
results will be given only for A = H = 0. 

In order now to satisfy Eq. as we have already 
stated, we impose equality of the last terms of Eqs. ®. 
This can be accomplished if we choose the following con¬ 
ditions to be satisfied 


P^MV^^ = PgUiV^^ 


Vg^ = 


IMP, 

mPe 




I mPe 


V 

Rl > 2p,MV,lRl 


a. 


( 12 ) 


where Rmax is the maximum value of R, to be specified 
below. We set P := 2piMV^^ = 2pgmV^^. 


Using fi, fg in Eqs. (g]) and then Eq. ® we obtain 

exp{-e{Pi + Pg)<^) = E{R, T) := ^ 

-ti 


Eg := 1 


F, := 1 


P!v^%e^^yRl 
v^i - PiRyElf""^ [ 1 - PiRVRl) 


Ck-i 


v^i - PIRPIRI)""^ \ 1 - P{RVRl) 


(13) 

Thus we verify that in general $ = <I>(i?, dt). From this 
we find that 

fii 

exp{—ePi^) = F 

exp{ePg^) = F (14) 

Computing the current density of Eq. @ we obtain 

r . ‘^mleR'f/Rl 

J(f) = enoaiexp[—epi^)— -x 


X exp 


[l-P{RyRl)]F^ 

'pfv^y^yRl' 


+ 


+ enoagexp{ePg^) 


1 - PiR^Rl) ^ 

2PgV^%eR^lRl 


X exp 


[1 - /3(i?Vit!2)]3/2 
l3^eVg%e^^yRl' 

1 - P(RPIRI 


(15) 


where Eqs. m and (HI are used. 

The magnetic field can be written as 

B = |e^ + V4/x^ 


(16) 


where lo/R is the magnitude of a vacuum toroidal field 
at some value of R. Projecting the curl of B on and 
equating it with the current density we obtain the GS- 
like equation to be specified and solved numerically in 
Section III. 


III. THE GRAD-SHAFRANOV-LIKE 
EQUATION 


For ITER-like equilibria we have Rq = 6.2m for the 
major radius and a = 2m the minor radius and therefore 
the aspect ratio is eg = 0.32. Thus the maximum distance 
perpendicular to the axis of symmetry is found to be 
Rmax = 7?o(l + eo)- Introducing the normalized variables 
p := R/Ro, C ■= z/Rq they range as 0.7 < p < 1.2 
and Id < 0.6. Assigning to the temperatures the values 
kTs ~ lOkeV, {s = i,e), we find P ~ 0.021. Thus we 
must have Rn > Rmax'/P — 1.2m. We choose Rn = 2.m 
and define p„ := Rn/Ro- Also we take as typical values 
no — 10^®m“^ and — lO^m/sec. 
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From the restriction 0 < < 1 adopted it fol¬ 

lows Oe ~ 0.023ai. We proceed now to a normaliza¬ 
tion of all the above involved quantities. In particular, 
the poloidal flux is normalized as 'i>n ■= where 

dto := Rn/iPiVi^e) ~ 0.2VF5. Then we have 

^ (l-/3(PVPn)) 

with the form of Eqs. (HI remaining unaffected. Subse¬ 
quently we obtain 


J0 = [Ciaiexp{-el}i^) 

HipVpD-^r. 


CetteeXp^ePe^)] X 

Vl/2 




(18) 


The numerical values of the constants are Ci ~ 7405 
and Ce = {M/m)Ci. The completely normalized GS-like 
equation turns out to be 


dp"^ p dp dC,'^ 


(19) 


where $(i?, 4>„) is determined by Eq. (|T^ . It is noted 
that (IT^ holds for arbitrary values of A and B. 


IV. EQUILIBRIA AND EQUILIBRIUM 
PROPERTIES 


First we define the fixed boundary coinciding with out¬ 
ermost magnetic surface as follows. The equation for the 
upper part of the bounding surface, which if taken to 
hold for the lower part as well would give a symmetric 
boundary, is 


Pb = 1 + eocos{T + asin{T)) 

Cb — CmaxSlTliT^ (^ 0 ) 


where Cmax = with 5 = (1—p 5 )/eo, and a = sin~^(S). 
Thus the following relations hold: ps = 1 — Seo and 0s = 
TT — tan~^{K/S) (see Fig. 1) The parameter r is any 
increasing function of the polar angle 0, satisfying r(0) = 
0, t(tt) = TT and t{9s) = tt/2. In our model we take 

t{ 9) = to9^+ti9^ 


7r0" - 027r"-i 
-01 + 

7r0^ - 0^7r"-i 


( 21 ) 


with n = 8. In order to complete the asymmetric bound¬ 
ing curve we specify now the lower part of it (C < 0) as 
follows. The left lower branch of the curve is given by 


Pf, = 1 -I- eocos{9) 


Cb 

Pi 


-[2pieo(l -I- COS0)] 


2eo(l -I- cos9s) ’ 


1/2 

[tt < 9 < 2n — 9s) 


( 22 ) 


while the right lower branch of the curve is given by 


Pb 

Cb 

P2 


1 -I- eocos(0) 

-[2p2eo(l - cos0)]^/^ 

_ Smax _ f 27 r 

2eo(l — cos9s )’ 


9s <9 < 2 tt) (23) 


The divertor X-point for the asymmetric equilibrium is 
located at px = l + eocos9s = 0.9139 and Cx = —Cmax = 
-0.6105. 



FIG. 1: The boundary determined by the parametric Eqs. 
(20)-(23) 

We adopt the numerical integration scheme described 
in detail in Sec. 4 of our previous paper [l^ using the 
nine point formula. Using the numerical algorithm asso¬ 
ciated with this formula for peaked toroidal current den¬ 
sity (A = B = —1) we obtained the equilibrium shown 
in Fig. 2. The magnetic axis was found by the numerical 
procedure to be located at (paXa) = (1.0574,0.02). The 
numerical procedure converged to the actual solution of 
Fig. 2 after N = 2027 iterations with a finite differ¬ 
ence step size of /i = 0.01. We have taken ai = 0.0015. 
At the boundary we have set 4>„(6) = 1.0, while the 
numerical integration scheme turned out to give for the 
magnetic axis 4'„(a) = 0.0. The configuration of Fig. 
2 remains nearly unaffected when changes to hollow 
{A = B = 0). 

On the basis of the solution constructed we examined 
the characteristics for both equilibria with peaked and 
hollow by calculating the safety factor and certain 
equilibrium quantities as follows. For axisymmetric plas¬ 
mas the safety factor can be put in the form 

1 

9('Pn) = ^y^ Q{'^n,9)d9 (24) 
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FIG. 2: The equilibrium configuration in connection with the 
numerical solution of the GS-like equation [Eqs. (17)-(18)] 
for peaked toroidal current density A = B = —1 [Eq. (Illll l. 
The magnetic surfaces are nearly the same for the respective 
equilibrium with hollow (A = B = 0). 


where 


Q{^n,d) 


Qo 


[jp-ir+e] 

\{P - l)('I'ra).p + C('f'ra),cl 


(25) 


and Qq ~ 2.28 is a dimensionless constant. The variation 
of the safety factor for peaked increasing monotoni- 
cally from the magnetic axis to the plasma boundary is 
shown in Fig. 3. The profile of on the midplane 
C = 0 which vanishes on the boundary is shown in Fig. 
4(above), while the respective hollow, current-hole like 
profile of is shown in Fig. 4(below). Since in the later 
case vanishes on axis the safety factor tends to infin¬ 
ity thereon and the equilibrium in the central region has 
negative magnetic shear. Current hole equilibria of this 
kind have been observed in JET [T^ and JT-60U (l^ . 

The toroidal ion fluid velocity Ui^ is given by 




v^fcfv 


V7^adp/Pn)'&n ( K \ 
[i-/3(pVp^)]^^=^ ^ Vi-d(pVp ^)) 


1 + 


•y/l-/3(p2/p2) 




(26) 


This for the two equilibria considered is plotted in Fig. 5. 
Note that the profile shapes are similar to the respective 
shapes of the current density profiles. Associated electric 
field profiles are given in Fig. 6. 

According to 0 (p. 14 therein) the total pressure 


FIG. 3: Variation of the safety factor for the equilibrium of 
Fig. 2 from the magnetic axis to the boundary for peaked on 
axis J^. At the magnetic axis we have qa = 1.17. 


tensor can be defined as 


Pki = Pli+P^i = 

= M J{vk- U^k){vi - Uii)(fv + 

+ rn {vk- Uek)(vi - Uei)d^v, (fc, I = r,(j), z) 

(27) 

where P^i are the pressure tensors of the ion and 
electron fluids and Ui^, Ue^, are the ion and electron fluid 
velocities given by Ui^ = (l/n^) J v^fiCpv and Ue^ = 
(l/rie) / v^feCPv. In the limit of — > 0 we recover 

the usual formula of nok(Ti + Tg) for the total pressure. 
For the distribution functions (US the tensor becomes 
diagonal and anisotropic and we have Prr = Pzz '■= Pi 
and Prj)rj) := P 2 - Here 


Pi := Mli + mie = M 


vlfid^v + m 


/ 


vlfed^V 


P2 


MJi + mJe = M y (u0 - Ui^Yfid^v + 

m hv^- Ue 4 >ffed^V 

(28) 


Computing the integrals A, A, Ji, Je we calcualted Pi, 
P 2 and the index of anisotropy 


Lo := 



(29) 


The calculation though lengthy and tedious is straight¬ 
forward. According to the results shown in Fig. 7, Pi 



















{,Lu/v>i) r {pUj/v>i) 
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FIG. 4: Peaked on axis toroidal current density for A = B = 
— 1 in connection with the equilibrium of Fig. 2. (above) and 
respective hollow for A — B = 0 (below). 

and P 2 decrease very weakly from the magnetic axis to 
the plasma boundary, irrespective of the shape of the cur¬ 
rent density profile.Similar is the behavior of the density 
variation given in Fig. 8 but the anisotropy index shown 
in Fig. 9 is different; the latter follows the respective 
variation of the current density profile. Apparently, such 
nearly flat pressure and density profiles are not repre¬ 
sentative for tokamaks. They just may be regarded as 
an approximation of the nearly flat pressure and density 
profiles observed during the L-FI transition in the cen¬ 
tral region inside the edge pedestal at which the pressure 
drops sharply giving rise to the ELMs activity. This edge 
region can not be described by the equilibrium solutions 




FIG. 5: The ion fluid velocity profile on the midplane z = 0 
for the equilibrium with peaked (above) and hollow 
(below). 


constructed here; in particular, we had a difficulty to 
make Pi and P 2 vanish on the boundary. This in addi¬ 
tion to numerical reasons should be related with the fact 
that Vlasov theory involves particle orbits which, unlike 
fluid theories, are not compatible with a fixed boundary; 
in this sense Vlasov is not appropriate to describe fixed 
boundary equilibria in the region close to the boundary. 









|E_| (kV/m) 
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P(?=0) 


FIG. 6 : The electric field profile on the midplane z = 0 for the 

equilibrium with peaked (above) and hollow (below). 7 . pressures Pi and P 2 for the equilibrium 

with peaked J<^. Similar in shape are the respective pressure 
profiles with hollow J^. 

V. SUMMARY 


In the framework of Maxwell-Vlasov theory we have 
derived a transcendental GS-like equation [Eqs. (17)- 
(18)] for quasineutral axisymmetric plasmas by select¬ 
ing exponentially deformed Maxwellian ion and electron 
distribution functions depending on the two known con¬ 
stants of particle motion [Eqs. (10)]. Then we have 
solved this equation numerically for a tokamak pertinent 
fixed diverted plasma boundary. To avoid unphysical 
perpendicular drifts on axis the electric filed was pro¬ 
visioned to vanish thereon. Depending on the symme¬ 
try properties of the distribution functions in connec¬ 
tion with the values of pertinent free parameters, we 
derived equlibria with toroidal current density profiles 


either peaked on axis or hollow. Both equilibria have 
nearly the same magnetic surfaces, sheared toroidal ion 
flow and diagonal anisotropic pressure tensor with dif¬ 
ferent toroidal and poloidal elements. The profiles of 
the toroidal velocity and of the pressure isotropy index 
are similar to the respective profiles of the current den¬ 
sity. However, the profiles of the pressure elements being 
nearly flat just may approximately represent respective 
experimental profiles in the central region developed dur¬ 
ing the L-H transition. 

It is interesting to pursue obtaining other Vlasov equi¬ 
libria with alternative choices of distribution functions, 
in particular distribution functions potentially creating 









mass density (x 1.67 • 10 ° Kg/m'^) 



0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 


P (C=0) 

FIG. 8: Density profile for the equilibrium with peaked J^. 
Similar in shape is the density profile for the equilibrium with 
hollow J^. 

more realistic pressure profiles. This may require nu¬ 
merical integrations in the velocity space which coupled 
with spatial integrations constitute a challenging prob¬ 
lem. Also, since poloidal flows play a role in the tran¬ 
sitions to improved confinement regimes in tokamaks it 
is desirable that these equilibria involve flows of arbi¬ 
trary direction. However, as already mentioned in Sec. 
I, this remains a tough open problem requiring additional 
Vlasov constants of motion. 
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FIG. 9: The pressure anisotropy index for the equilibrium 
with peaked (above) and hollow (below). 
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